package scu.maqiang.numeric;

import static java.lang.Math.*;

/**
 * 计算几何类，提供在有限元计算过程中处理几何形状构型相关计算
 * 
 * @author 马强
 *
 */
public class CG {

	/**
	 * 计算二维与三维点与线段之间的距离
	 * @param P
	 * @param AB
	 * @return
	 */
	public static double distancePointToLine(double[] P, double[][] AB) {
		double[] A = AB[0];
		double[] B = AB[1];
		double lengthAB = MVO.distance(A, B);
		double[] vecAB = MVO.add(B, -1.0, A);
		double[] vecAP = MVO.add(P, -1.0, A);
		double alpha = MVO.dot_product(vecAB, vecAP) / lengthAB;
		if (alpha < 0) {
			return MVO.distance(P, A);
		}
		if (alpha > lengthAB) {
			return MVO.distance(P, B);
		}
		double[] PC = MVO.add(vecAP, -alpha / lengthAB, vecAB);
		return MVO.L2Norm(PC);
	}

	public static double distancePointToTriangle(double[] P, double[][] ABC) {
		double[] A = ABC[0];
		double[] B = ABC[1];
		double[] C = ABC[2];
		double[] vecAB = MVO.add(B, -1.0, A);
		double[] vecAC = MVO.add(C, -1.0, A);
		double[] vecAP = MVO.add(P, -1.0, A);
		double[] v = MVO.cross_product(vecAB, vecAC);
		double lenV = MVO.L2Norm(v);
		double alpha = MVO.dot_product(vecAP, v) / lenV;
		double[] N = new double[3];
		T3ShapeFunction(ABC, P,  N);
		double a = N[0];
		double b = N[1];
		double c = N[2];

		if (a < 0) {
			if (b < 0) {
				double[] vecCP = MVO.add(P, -1.0, C);
				return MVO.L2Norm(vecCP);
			} else if (c < 0) {
				double[] vecBP = MVO.add(P, -1.0, B);
				return MVO.L2Norm((vecBP));
			} else {
				double[][] lineBC = {B, C};
				return distancePointToLine(P, lineBC);
			}
		}
		if(b < 0) {
			if (c < 0) {
				return MVO.L2Norm(vecAP);
			} else {
				double[][] lineCA = {C, A};
				return distancePointToLine(P, lineCA);
			}
		}
		if (c < 0) {
			double[][] lineAB = {A, B};
			return distancePointToLine(P, lineAB);
		}
		return alpha;
	}

	public static double distanceLineToLine(double[][] AB, double[][] CD) {
		return 0.0;
	}


	public static double length2DL2(double[][] xy) {
		return Math.sqrt((xy[0][0] - xy[1][0]) * (xy[0][0] - xy[1][0]) + (xy[0][1] - xy[1][1]) * (xy[0][1] - xy[1][1]));
	}

	/**
	 * 给定三角形三顶点坐标，计算二维三角形面积
	 * 
	 * @param xy 3x2 数组，行i表示第i个顶点，列i表示该顶点第j个坐标值
	 * @return 三角形面积，如果顶点坐标逆时针排列，面积为正; 如果顺时针排列，面积为负
	 */
	public static double area2DT3(double[][] xy) {
		double[] x = { xy[0][0], xy[1][0], xy[2][0] };
		double[] y = { xy[0][1], xy[1][1], xy[2][1] };
		return 0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]));
	}

	public static Complex area2DT3(Complex[][] xy) {
		Complex[] x = { xy[0][0], xy[1][0], xy[2][0] };
		Complex[] y = { xy[0][1], xy[1][1], xy[2][1] };
		return ((x[1].sub(x[0])).mul(y[2].sub(y[0])).sub((x[2].sub(x[0])).mul(y[1].sub(y[0])))).mul(0.5);
	}

	public static int area2DT3(int[][] xy) {
		int[] x = { xy[0][0], xy[1][0], xy[2][0] };
		int[] y = { xy[0][1], xy[1][1], xy[2][1] };
		return ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0])) / 2;
	}


	
	public static int area2DT3(int[] p1, int[] p2, int[] p3) {
		return ((p2[0] - p1[0]) * (p3[1] - p1[1]) - (p3[0] - p1[0]) * (p2[1] - p1[1])) / 2;
	}
	
	public static double area2DT3(double[] p1, double[] p2, double[] p3) {
		return 0.5 * ((p2[0] - p1[0]) * (p3[1] - p1[1]) - (p3[0] - p1[0]) * (p2[1] - p1[1]));
	}
	/**
	 * 使用海伦公式计算三维三角形面积
	 * 
	 * @param xyz
	 *            3x3ά���飬��ά���������������, ��ָ���ʾ����ţ���ָ���ʾ������
	 * @return ���������
	 */
	public static double area3DT3(double[][] xyz) {
		double[] x = { xyz[0][0], xyz[1][0], xyz[2][0] };
		double[] y = { xyz[0][1], xyz[1][1], xyz[2][1] };
		double[] z = { xyz[0][2], xyz[1][2], xyz[2][2] };
		double a = sqrt((x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]) + (z[0] - z[1]) * (z[0] - z[1]));
		double b = sqrt((x[1] - x[2]) * (x[1] - x[2]) + (y[1] - y[2]) * (y[1] - y[2]) + (z[1] - z[2]) * (z[1] - z[2]));
		double c = sqrt((x[2] - x[0]) * (x[2] - x[0]) + (y[2] - y[0]) * (y[2] - y[0]) + (z[2] - z[0]) * (z[2] - z[0]));
		double p = 0.5 * (a + b + c);
		return sqrt(p * (p - a) * (p - b) * (p - c));
	}

	/**
	 * 计算二维四边形面积
	 * 
	 * @param xy
	 *            4x2 数组, 四边形结点坐标
	 * @return 面积
	 */
	public static double area2DQ4(double[][] xy) {
		double[][] xy1 = { xy[0], xy[1], xy[2] };
		double[][] xy2 = { xy[0], xy[2], xy[3] };
		return area2DT3(xy1) + area2DT3(xy2);
	}

	/**
	 * 计算三维四边形面积
	 * @param xyz
	 * @return
	 */
	public static double area3DQ4(double[][] xyz) {
		double[][] xyz1 = { { xyz[0][0], xyz[1][0], xyz[2][0] }, { xyz[0][1], xyz[1][1], xyz[2][1] },
				{ xyz[0][2], xyz[1][2], xyz[2][2] } };
		double[][] xyz2 = { { xyz[0][0], xyz[2][0], xyz[3][0] }, { xyz[0][1], xyz[2][1], xyz[3][1] },
				{ xyz[0][2], xyz[2][2], xyz[3][2] } };
		return area3DT3(xyz1) + area3DT3(xyz2);
	}

	/**
	 * 计算三维四面体体积
	 * 
	 * @param coord 4x3实数组, 存储四面体四个顶点坐标
	 * @return 四面体体积
	 */
	public static double volume3DT4(double[][] coord) {
		double[] x = MVO.col(coord, 0);
		double[] y = MVO.col(coord, 1);
		double[] z = MVO.col(coord, 2);
		double x1x0 = x[1] - x[0];
		double x2x0 = x[2] - x[0];
		double x3x0 = x[3] - x[0];
		double y1y0 = y[1] - y[0];
		double y2y0 = y[2] - y[0];
		double y3y0 = y[3] - y[0];
		double z1z0 = z[1] - z[0];
		double z2z0 = z[2] - z[0];
		double z3z0 = z[3] - z[0];
		double volume6 = x1x0 * (y2y0 * z3z0 - y3y0 * z2z0) - x2x0 * (y1y0 * z3z0 - y3y0 * z1z0)
				+ x3x0 * (y1y0 * z2z0 - y2y0 * z1z0);
		return Math.abs(volume6) / 6.0;
	}

	public static double volume3DT4_Signed(double[][] coord) {
		double[] x = MVO.col(coord, 0);
		double[] y = MVO.col(coord, 1);
		double[] z = MVO.col(coord, 2);
		double x1x0 = x[1] - x[0];
		double x2x0 = x[2] - x[0];
		double x3x0 = x[3] - x[0];
		double y1y0 = y[1] - y[0];
		double y2y0 = y[2] - y[0];
		double y3y0 = y[3] - y[0];
		double z1z0 = z[1] - z[0];
		double z2z0 = z[2] - z[0];
		double z3z0 = z[3] - z[0];
		double volume6 = x1x0 * (y2y0 * z3z0 - y3y0 * z2z0) - x2x0 * (y1y0 * z3z0 - y3y0 * z1z0)
				+ x3x0 * (y1y0 * z2z0 - y2y0 * z1z0);
		return volume6 / 6.0;
	}

	/**
	 * 计算三维六面体体积
	 * 
	 * @param xyz
	 *            8x3实数组, 存储六面体8个顶点坐标
	 * @return 六面体体积
	 */
	public static double volume3DH8(double[][] xyz) {
		double[][] xyz1 = { xyz[0], xyz[1], xyz[2], xyz[6] };
		double result = volume3DT4(xyz1);
		double[][] xyz2 = { xyz[0], xyz[5], xyz[1], xyz[6] };
		result += volume3DT4(xyz2);
		double[][] xyz3 = { xyz[0], xyz[2], xyz[3], xyz[6] };
		result += volume3DT4(xyz3);
		double[][] xyz4 = { xyz[0], xyz[4], xyz[5], xyz[6] };
		result += volume3DT4(xyz4);
		double[][] xyz5 = { xyz[0], xyz[3], xyz[7], xyz[6] };
		result += volume3DT4(xyz5);
		double[][] xyz6 = { xyz[0], xyz[7], xyz[4], xyz[6] };
		result += volume3DT4(xyz6);
		return result;
	}

	/**
	 * 计算三角形三个角的角度值
	 * @param xy
	 * @param angle
	 */
	public static void angle2DT3(double[][] xy, double[] angle) {
		double[] x = { xy[0][0], xy[1][0], xy[2][0] };
		double[] y = { xy[0][1], xy[1][1], xy[2][1] };
		double[] length = new double[3];
		length[0] = Math.sqrt((x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]));
		length[1] = Math.sqrt((x[1] - x[2]) * (x[1] - x[2]) + (y[1] - y[2]) * (y[1] - y[2]));
		length[2] = Math.sqrt((x[0] - x[2]) * (x[0] - x[2]) + (y[0] - y[2]) * (y[0] - y[2]));

		angle[0] = Math.acos(((x[1] - x[0]) * (x[2] - x[0]) + (y[1] - y[0]) * (y[2] - y[0])) / (length[0] * length[2]));
		angle[1] = Math.acos(((x[2] - x[1]) * (x[0] - x[1]) + (y[2] - y[1]) * (y[0] - y[1])) / (length[0] * length[1]));
		angle[2] = Math.acos(((x[0] - x[2]) * (x[1] - x[2]) + (y[0] - y[2]) * (y[1] - y[2])) / (length[1] * length[2]));
	}

	public static void angle2DT6(double[][] xy, double[] angle) {
		angle2DT3(xy, angle);
		angle[3] = Math.acos(-1.0);
		angle[4] = angle[3];
		angle[5] = angle[3];
	}

	public static void angle2DQ4(double[][] xy, double[] angle) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] length = new double[4];
		length[0] = Math.sqrt((x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]));
		length[1] = Math.sqrt((x[1] - x[2]) * (x[1] - x[2]) + (y[1] - y[2]) * (y[1] - y[2]));
		length[2] = Math.sqrt((x[2] - x[3]) * (x[2] - x[3]) + (y[2] - y[3]) * (y[2] - y[3]));
		length[3] = Math.sqrt((x[3] - x[0]) * (x[3] - x[0]) + (y[3] - y[0]) * (y[3] - y[0]));

		angle[0] = Math.acos(((x[1] - x[0]) * (x[3] - x[0]) + (y[1] - y[0]) * (y[3] - y[0])) / (length[0] * length[3]));
		angle[1] = Math.acos(((x[2] - x[1]) * (x[0] - x[1]) + (y[2] - y[1]) * (y[0] - y[1])) / (length[1] * length[0]));
		angle[2] = Math.acos(((x[3] - x[2]) * (x[1] - x[2]) + (y[3] - y[2]) * (y[1] - y[2])) / (length[2] * length[1]));
		angle[3] = Math.acos(((x[0] - x[3]) * (x[2] - x[3]) + (y[0] - y[3]) * (y[2] - y[3])) / (length[3] * length[2]));
	}

	public static void angle2DQ8(double[][] xy, double[] angle) {
		angle2DQ4(xy, angle);
		angle[4] = angle[5] = angle[6] = angle[7] = Math.acos(-1.0);
	}

	public static void edge2DT3(double[][] xy, double[] length) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] dx = {x[0] - x[1], x[1] - x[2], x[2] - x[0]};
		double[] dy = {y[0] - y[1], y[1] - y[2], y[2] - y[0]};
		for(int i = 0; i < 3; i++) {
			length[i] = Math.sqrt(dx[i] * dx[i] + dy[i] * dy[i]);
		}
	}

	public static void edge2DT3(double[][] xy, double[] length, double[] nx, double[] ny) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] dx = {x[0] - x[1], x[1] - x[2], x[2] - x[0]};
		double[] dy = {y[0] - y[1], y[1] - y[2], y[2] - y[0]};
		
		for(int i = 0; i < 3; i++) {
			length[i] = Math.sqrt(dx[i] * dx[i] + dy[i] * dy[i]);
			nx[i] = dy[i] / length[i];
			ny[i] = -dx[i] / length[i];
		}
	}
	
	public static void length2DT3(double[][] xy, double[] length) {
		double[] x = { xy[0][0], xy[1][0], xy[2][0] };
		double[] y = { xy[0][1], xy[1][1], xy[2][1] };
		length[0] = Math.sqrt((x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]));
		length[1] = Math.sqrt((x[1] - x[2]) * (x[1] - x[2]) + (y[1] - y[2]) * (y[1] - y[2]));
		length[2] = Math.sqrt((x[0] - x[2]) * (x[0] - x[2]) + (y[0] - y[2]) * (y[0] - y[2]));
	}
	
	public static boolean pointIn1DL2(double[][] x, double[] x0, double[] N) {
		double[] xx = MVO.col(x, 0);
		double xx0 = x0[0];
		double length = Math.abs(xx[0] - xx[1]);
		double w1 = Math.abs(xx[1] - xx0) / length;
		double w2 = Math.abs(xx0 - xx[0]) / length;
		if (Math.abs(1.0 - w1 - w2) < Constants.Er) {
			N[0] = w1;
			N[1] = w2;
			return true;
		} else {
			return false;
		}
	}
	
	public static boolean pointInT3(double[][] xy, double[] xy0, double[] N) {
		CG.T3ShapeFunction(xy, xy0, N);
		return Math.abs(MVO.sum(N) - 1.0) < 1.0e-10;
	}

	public static boolean pointInT4(double[][] xy, double[] xy0, double[] N) {
		CG.T4ShapeFunction(xy, xy0, N);
		return Math.abs(MVO.sum(N) - 1.0) < 1.0e-10;
	}
	
	public static boolean pointInQ4(double[][] xy, double[] xy0, double[] N) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double area0 = area2DQ4(xy);

		double[][] newXY = new double[3][];
		newXY[0] = new double[] { xy0[0], xy0[1] };
		newXY[1] = new double[] { x[0], y[0] };
		newXY[2] = new double[] { x[1], y[1] };
		double area1 = Math.abs(area2DT3(newXY));

		newXY[1] = new double[] { x[1], y[1] };
		newXY[2] = new double[] { x[2], y[2] };
		double area2 = Math.abs(area2DT3(newXY));

		newXY[1] = new double[] { x[2], y[2] };
		newXY[2] = new double[] { x[3], y[3] };
		double area3 = Math.abs(area2DT3(newXY));

		newXY[1] = new double[] { x[3], y[3] };
		newXY[2] = new double[] { x[0], y[0] };
		double area4 = Math.abs(area2DT3(newXY));
		
		if(Math.abs(area0 - area1 - area2 - area3 - area4) < 1.0e-10) {
			CG.Q4ShapeFunction(xy, xy0, N);
			return true;
		}
		return false;
	}
	
	public static boolean pointInQ4(double[][] xy, double[] xy0) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double area0 = area2DQ4(xy);

		double[][] newXY = new double[3][];
		newXY[0] = new double[] { xy0[0], xy0[1] };
		newXY[1] = new double[] { x[0], y[0] };
		newXY[2] = new double[] { x[1], y[1] };
		double area1 = Math.abs(area2DT3(newXY));

		newXY[1] = new double[] { x[1], y[1] };
		newXY[2] = new double[] { x[2], y[2] };
		double area2 = Math.abs(area2DT3(newXY));

		newXY[1] = new double[] { x[2], y[2] };
		newXY[2] = new double[] { x[3], y[3] };
		double area3 = Math.abs(area2DT3(newXY));

		newXY[1] = new double[] { x[3], y[3] };
		newXY[2] = new double[] { x[0], y[0] };
		double area4 = Math.abs(area2DT3(newXY));
		
		return Math.abs(area0 - area1 - area2 - area3 - area4) < 1.0e-10;
	}
	
	public static boolean pointInSquare(double[][] xy, double[] xy0) {
		double[] xyL = xy[0];
		double[] xyR = xy[2];
		return (xy0[0]>= xyL[0] && xy0[0] <= xyR[0]) && (xy0[1]>= xyL[1] && xy0[1] <= xyR[1]);
	}

	public static boolean pointInH8(double[][] xyz, double[] xyz0) {
		double volume0 = volume3DH8(xyz);

		double volume = 0.0;
		double[][] newXYZ1 = { xyz[0], xyz[1], xyz[3], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ1));

		double[][] newXYZ2 = { xyz[1], xyz[2], xyz[3], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ2));

		double[][] newXYZ3 = { xyz[4], xyz[7], xyz[5], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ3));

		double[][] newXYZ4 = { xyz[7], xyz[6], xyz[5], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ4));

		double[][] newXYZ5 = { xyz[0], xyz[3], xyz[7], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ5));

		double[][] newXYZ6 = { xyz[0], xyz[7], xyz[4], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ6));

		double[][] newXYZ7 = { xyz[1], xyz[5], xyz[6], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ7));

		double[][] newXYZ8 = { xyz[1], xyz[6], xyz[2], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ8));

		double[][] newXYZ9 = { xyz[0], xyz[5], xyz[1], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ9));

		double[][] newXYZ10 = { xyz[0], xyz[4], xyz[5], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ10));

		double[][] newXYZ11 = { xyz[2], xyz[6], xyz[3], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ11));

		double[][] newXYZ12 = { xyz[3], xyz[6], xyz[7], xyz0 };
		volume += Math.abs(volume3DT4(newXYZ12));

		return Math.abs(volume0 - volume) < 1.0e-10;
	}

	public static boolean pointInCube(double[][] xyz, double[] xyz0) {
		double[] xyL = xyz[0];
		double[] xyR = xyz[6];
		return (xyz0[0]>= xyL[0] && xyz0[0] <= xyR[0]) && (xyz0[1]>= xyL[1] && xyz0[1] <= xyR[1]) && (xyz0[2]>= xyL[2] && xyz0[2] <= xyR[2]);
	}
	/**
	 * ������ά�����������������
	 * 
	 * @param xyz
	 *            4x3���飬������4�����
	 * @param sAngle
	 *            4ά���飬������4������������
	 */
	public static void solidAngle3DT4(double[][] xyz, double[] sAngle) {
		// double[] x = MVO.col(xyz, 0);
		// double[] y = MVO.col(xyz, 1);
		// double[] z = MVO.col(xyz, 2);
	}

	public static double[][] deformGradient2DT3(double[] X, double[] Y, double[] x, double[] y) {
		double[][] dxds = {{-x[0] + x[1], -x[0] + x[1]}, {-y[0] + y[1], -y[0] + y[2]}};
		double[][] dXds = {{-X[0] + X[1], -X[0] + X[1]}, {-Y[0] + Y[1], -Y[0] + Y[2]}};
		double[][] dsdX = MVO.inverseMat2(dxds);
		return MVO.matmul(dxds, dsdX);
	}

	public static double gradient2DT3(double[] x, double[] y, double[][] grad) {
		double area = 0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]));
		double A2D = 1 / (2 * area);
		grad[0][0] = A2D * (y[1] - y[2]);
		grad[1][0] = A2D * (y[2] - y[0]);
		grad[2][0] = A2D * (y[0] - y[1]);
		grad[0][1] = A2D * (x[2] - x[1]);
		grad[1][1] = A2D * (x[0] - x[2]);
		grad[2][1] = A2D * (x[1] - x[0]);
		return area;
	}

	/**
	 * 计算二维三角形单元形函数梯度
	 * 
	 * @param coord
	 *            3x2 数组, 存储单元坐标
	 * @param grad
	 *            3x2ά���飬�����θ�����κ�����x��y�����ݶ�ֵ
	 * @return ���������
	 */
	public static double gradient2DT3(double[][] coord, double[][] grad) {
		double[] x = { coord[0][0], coord[1][0], coord[2][0] };
		double[] y = { coord[0][1], coord[1][1], coord[2][1] };
		double area = 0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]));
		double A2D = 1 / (2 * area);
		grad[0][0] = A2D * (y[1] - y[2]);
		grad[1][0] = A2D * (y[2] - y[0]);
		grad[2][0] = A2D * (y[0] - y[1]);
		grad[0][1] = A2D * (x[2] - x[1]);
		grad[1][1] = A2D * (x[0] - x[2]);
		grad[2][1] = A2D * (x[1] - x[0]);
		return area;
	}

	public static Complex gradient2DT3(Complex[][] coord, Complex[][] grad) {
		Complex[] x = { coord[0][0], coord[1][0], coord[2][0] };
		Complex[] y = { coord[0][1], coord[1][1], coord[2][1] };
		Complex area = (((x[1].sub(x[0])).mul(y[2].sub(y[0])).sub(x[2].sub(x[0]).mul(y[1].sub(y[0]))))).mul(0.5);
		Complex A2D = area.mul(2.0);
		grad[0][0] = (y[1].sub(y[2])).div(A2D);
		grad[1][0] = (y[2].sub(y[0])).div(A2D);
		grad[2][0] = (y[0].sub(y[1])).div(A2D);
		grad[0][1] = (x[2].sub(x[1])).div(A2D);
		grad[1][1] = (x[0].sub(x[2])).div(A2D);
		grad[2][1] = (x[1].sub(x[0])).div(A2D);
		return area;
	}


	// public static double gradient2DQ4(double[][] coord, double[][] grad) {
	// double[][] coef = new double[4][4];
	// Q4ShapeFuncionCoef(coord, coef);
	// double area = area2DQ4(coord);
	// for(int i = 0; i < 4; i++) {
	// grad[i][0] = coef[i][1] + coef[i][3] * coord[i][1];
	// grad[i][1] = coef[i][2] + coef[i][3] * coord[i][0];
	// }
	// return area;
	// }

	public static double gradient3DT4(double[][] coord, double[][] grad) {
		double[] x = MVO.col(coord, 0);
		double[] y = MVO.col(coord, 1);
		double[] z = MVO.col(coord, 2);
		double x1x0 = x[1] - x[0];
		double x2x0 = x[2] - x[0];
		double x3x0 = x[3] - x[0];
		double x2x1 = x[2] - x[1];
		double x3x1 = x[3] - x[1];
		double y1y0 = y[1] - y[0];
		double y2y0 = y[2] - y[0];
		double y3y0 = y[3] - y[0];
		double y2y1 = y[2] - y[1];
		double y3y1 = y[3] - y[1];
		double z1z0 = z[1] - z[0];
		double z2z0 = z[2] - z[0];
		double z3z0 = z[3] - z[0];
		double z2z1 = z[2] - z[1];
		double z3z1 = z[3] - z[1];

		double volume6 = x1x0 * (y2y0 * z3z0 - y3y0 * z2z0) - x2x0 * (y1y0 * z3z0 - y3y0 * z1z0)
				+ x3x0 * (y1y0 * z2z0 - y2y0 * z1z0);
		grad[0][0] = -(y2y1 * z3z1 - y3y1 * z2z1) / volume6;
		grad[0][1] = +(x2x1 * z3z1 - x3x1 * z2z1) / volume6;
		grad[0][2] = -(x2x1 * y3y1 - x3x1 * y2y1) / volume6;
		grad[1][0] = +(y2y0 * z3z0 - y3y0 * z2z0) / volume6;
		grad[1][1] = -(x2x0 * z3z0 - x3x0 * z2z0) / volume6;
		grad[1][2] = +(x2x0 * y3y0 - x3x0 * y2y0) / volume6;
		grad[2][0] = -(y1y0 * z3z0 - y3y0 * z1z0) / volume6;
		grad[2][1] = +(x1x0 * z3z0 - x3x0 * z1z0) / volume6;
		grad[2][2] = -(x1x0 * y3y0 - x3x0 * y1y0) / volume6;
		grad[3][0] = +(y1y0 * z2z0 - y2y0 * z1z0) / volume6;
		grad[3][1] = -(x1x0 * z2z0 - x2x0 * z1z0) / volume6;
		grad[3][2] = +(x1x0 * y2y0 - x2x0 * y1y0) / volume6;

		return Math.abs(volume6 / 6);
	}

//	public static double ISOMap1D(ISOMapFunction1D isoFunc, double[] x, double xi, double[] v, double[] dx) {
//		int n = x.length;
//		double[] dXi = new double[n];
//		isoFunc.func(xi, v, dXi);
//		double j11 = MVO.dot_product(dXi, x);
//		double j12 = MVO.dot_product(dXi, y);
//		double j21 = MVO.dot_product(dEt, x);
//		double j22 = MVO.dot_product(dEt, y);
//		double detJ = j11 * j22 - j12 * j21;
//		for (int i = 0; i < n; i++) {
//			dx[i] = (j22 * dXi[i] - j12 * dEt[i]) / detJ;
//			dy[i] = (-j21 * dXi[i] + j11 * dEt[i]) / detJ;
//		}
//		return detJ;
//	}
	public static double ISOMap1D(ISOMapFunc1D isoFunc, double[] x, double xi, double[] v, double[] dx) {
		int n = x.length;
		double[] dXi = new double[n];
		isoFunc.func(xi, v, dXi);
		double detJ = MVO.dot_product(x, dXi);
		for(int i = 0; i < n; i++) {
			dx[i] = dXi[i] / detJ;
		}
		return detJ;
	}

	public static double ISOMap2D(ISOMapFunc2D isoFunc, double[] x, double[] y, double xi, double et, double[] v,
			double[] dx, double[] dy) {
		int n = x.length;
		double[] dXi = new double[n];
		double[] dEt = new double[n];
		isoFunc.func(xi, et, v, dXi, dEt);
		double j11 = MVO.dot_product(dXi, x);
		double j12 = MVO.dot_product(dXi, y);
		double j21 = MVO.dot_product(dEt, x);
		double j22 = MVO.dot_product(dEt, y);
		double detJ = j11 * j22 - j12 * j21;
		for (int i = 0; i < n; i++) {
			dx[i] = (j22 * dXi[i] - j12 * dEt[i]) / detJ;
			dy[i] = (-j21 * dXi[i] + j11 * dEt[i]) / detJ;
		}
		return detJ;
	}
	
	public static double ISOMap3D(ISOMapFunc3D isoFunc, double[] x, double[] y, double[] z, double xi, double et,
			double zt, double[] v, double[] dx, double[] dy, double[] dz) {
		int n = x.length;
		double[] dXi = new double[n];
		double[] dEt = new double[n];
		double[] dZt = new double[n];
		isoFunc.func(xi, et, zt, v, dXi, dEt, dZt);
		double j11 = MVO.dot_product(dXi, x);
		double j12 = MVO.dot_product(dXi, y);
		double j13 = MVO.dot_product(dXi, z);
		double j21 = MVO.dot_product(dEt, x);
		double j22 = MVO.dot_product(dEt, y);
		double j23 = MVO.dot_product(dEt, z);
		double j31 = MVO.dot_product(dZt, x);
		double j32 = MVO.dot_product(dZt, y);
		double j33 = MVO.dot_product(dZt, z);
		double detJ = j11 * (j22 * j33 - j32 * j23) - j12 * (j21 * j33 - j31 * j23) + j13 * (j21 * j32 - j31 * j22);
		double B11 = j22 * j33 - j32 * j23;
		double B21 = -(j21 * j33 - j31 * j23);
		double B31 = j21 * j32 - j31 * j22;
		double B12 = -(j12 * j33 - j32 * j13);
		double B22 = j11 * j33 - j31 * j13;
		double B32 = -(j11 * j32 - j31 * j12);
		double B13 = j12 * j23 - j22 * j13;
		double B23 = -(j11 * j23 - j21 * j13);
		double B33 = j11 * j22 - j21 * j12;
		for (int i = 0; i < n; i++) {
			dx[i] = (dXi[i] * B11 + dEt[i] * B12 + dZt[i] * B13) / detJ;
			dy[i] = (dXi[i] * B21 + dEt[i] * B22 + dZt[i] * B23) / detJ;
			dz[i] = (dXi[i] * B31 + dEt[i] * B32 + dZt[i] * B33) / detJ;
		}
		return detJ;
	}

	public static double ISOMap3DTo2D(ISOMapFunc2D isoFunc, double[] x, double[] y, double[] z, double xi, double et, double[] u) {
		int n = x.length;
		double[] dXi = new double[n];
		double[] dEt = new double[n];
		isoFunc.func(xi, et, u, dXi, dEt);
		double j11 = MVO.dot_product(dXi, x);
		double j12 = MVO.dot_product(dXi, y);
		double j13 = MVO.dot_product(dXi, z);
		double j21 = MVO.dot_product(dEt, x);
		double j22 = MVO.dot_product(dEt, y);
		double j23 = MVO.dot_product(dEt, z);
		double gx = j12 * j23 - j22 * j13;
		double gy = j21 * j13 - j11 * j23;
		double gz = j11 * j22 - j21 * j12;
		double detJ = sqrt(gx * gx + gy * gy + gz * gz);
		return detJ;
	}

	public static void T3ShapeFunction(double[][] xy, double[] xy0, double[] N) {
		double[] x = { xy[0][0], xy[1][0], xy[2][0] };
		double[] y = { xy[0][1], xy[1][1], xy[2][1] };

		double area0 = 0.5 * Math.abs((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]));
		double area1 = 0.5 * Math.abs((x[1] - xy0[0]) * (y[2] - xy0[1]) - (x[2] - xy0[0]) * (y[1] - xy0[1]));
		double area2 = 0.5 * Math.abs((xy0[0] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (xy0[1] - y[0]));
		double area3 = 0.5 * Math.abs((x[1] - x[0]) * (xy0[1] - y[0]) - (xy0[0] - x[0]) * (y[1] - y[0]));
		N[0] = area1 / area0;
		N[1] = area2 / area0;
		N[2] = area3 / area0;
	}

	public static double T3ShapeFunction(double[][] xy, double[] xy0, double[] N, double[][] gradN) {
		double[] x = { xy[0][0], xy[1][0], xy[2][0] };
		double[] y = { xy[0][1], xy[1][1], xy[2][1] };

		double area0 = 0.5 * Math.abs((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]));
		double area1 = 0.5 * Math.abs((x[1] - xy0[0]) * (y[2] - xy0[1]) - (x[2] - xy0[0]) * (y[1] - xy0[1]));
		double area2 = 0.5 * Math.abs((xy0[0] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (xy0[1] - y[0]));
		double area3 = 0.5 * Math.abs((x[1] - x[0]) * (xy0[1] - y[0]) - (xy0[0] - x[0]) * (y[1] - y[0]));
		N[0] = area1 / area0;
		N[1] = area2 / area0;
		N[2] = area3 / area0;
		double A2D = 1 / (2 * area0);
		gradN[0][0] = A2D * (y[1] - y[2]);
		gradN[1][0] = A2D * (y[2] - y[0]);
		gradN[2][0] = A2D * (y[0] - y[1]);
		gradN[0][1] = A2D * (x[2] - x[1]);
		gradN[1][1] = A2D * (x[0] - x[2]);
		gradN[2][1] = A2D * (x[1] - x[0]);
		return area0;
	}

	public static void T6ShapeFunction(double[][] xy, double[] xy0, double[] N) {
		double[] L = new double[3];
		T3ShapeFunction(xy, xy0, L);
		N[0] = L[0] * (2 * L[0] - 1);
		N[1] = L[1] * (2 * L[1] - 1);
		N[2] = L[2] * (2 * L[2] - 1);
		N[3] = 4 * L[0] * L[1];
		N[4] = 4 * L[1] * L[2];
		N[5] = 4 * L[0] * L[2];
	}

	/**
	 * 
	 * @param xy
	 *            6x2���飬6��������ε�Ԫ������
	 * @param L
	 *            3ά���飬��������������ĵ�Ԫ�κ���ֵ
	 * @param N
	 *            6ά���飬6����κ���ֵ
	 * @param gradN
	 *            6x2������6����κ�������ֵ
	 * @return
	 */
	public static double T6ShapeFunction(double[][] xy, double[] L, double[] N, double[][] gradN) {
		double[][] xxyy = MVO.subMatrix(xy, 0, 2, 0, 1);
		double[][] grad = new double[3][2];
		double area = gradient2DT3(xxyy, grad);
		double[] dxL = { grad[0][0], grad[1][0], grad[2][0] };
		double[] dyL = { grad[0][1], grad[1][1], grad[2][1] };

		double L1 = L[0], L2 = L[1], L3 = L[2];
		N[0] = L1 * (2 * L1 - 1);
		N[1] = L2 * (2 * L2 - 1);
		N[2] = L3 * (2 * L3 - 1);
		N[3] = 4 * L1 * L2;
		N[4] = 4 * L2 * L3;
		N[5] = 4 * L3 * L1;
		gradN[0][0] = 4 * dxL[0] * L1 - dxL[0];
		gradN[1][0] = 4 * dxL[1] * L2 - dxL[1];
		gradN[2][0] = 4 * dxL[2] * L3 - dxL[2];
		gradN[3][0] = 4 * (dxL[0] * L2 + L1 * dxL[1]);
		gradN[4][0] = 4 * (dxL[1] * L3 + L2 * dxL[2]);
		gradN[5][0] = 4 * (dxL[0] * L3 + L1 * dxL[2]);
		gradN[0][1] = 4 * dyL[0] * L1 - dyL[0];
		gradN[1][1] = 4 * dyL[1] * L2 - dyL[1];
		gradN[2][1] = 4 * dyL[2] * L3 - dyL[2];
		gradN[3][1] = 4 * (dyL[0] * L2 + L1 * dyL[1]);
		gradN[4][1] = 4 * (dyL[1] * L3 + L2 * dyL[2]);
		gradN[5][1] = 4 * (dyL[0] * L3 + L1 * dyL[2]);
		return area;
	}

	public static void T4ShapeFunction(double[][] xyz, double[] xyz0, double[] N) {

		double volume0 = Math.abs(volume3DT4(xyz));

		double[][] subxyz1 = { xyz[1], xyz[2], xyz[3], xyz0 };
		double volume1 = Math.abs(volume3DT4(subxyz1));

		double[][] subxyz2 = { xyz[0], xyz[2], xyz[3], xyz0 };
		double volume2 = Math.abs(volume3DT4(subxyz2));

		double[][] subxyz3 = { xyz[0], xyz[3], xyz[1], xyz0 };
		double volume3 = Math.abs(volume3DT4(subxyz3));

		double[][] subxyz4 = { xyz[0], xyz[1], xyz[2], xyz0 };
		double volume4 = Math.abs(volume3DT4(subxyz4));

		N[0] = volume1 / volume0;
		N[1] = volume2 / volume0;
		N[2] = volume3 / volume0;
		N[3] = volume4 / volume0;
	}

	public static void T10ShapeFunction(double[][] xyz, double[] xyz0, double[] N) {
		double[] L = new double[4];
		T4ShapeFunction(xyz, xyz0, L);
		double L1 = L[0], L2 = L[1], L3 = L[2], L4 = L[3];
		N[0] = L1 * (2 * L1 - 1);
		N[1] = L2 * (2 * L2 - 1);
		N[2] = L3 * (2 * L3 - 1);
		N[3] = L4 * (2 * L4 - 1);
		N[4] = 4 * L1 * L2;
		N[5] = 4 * L2 * L3;
		N[6] = 4 * L3 * L1;
		N[7] = 4 * L1 * L4;
		N[8] = 4 * L2 * L4;
		N[9] = 4 * L3 * L4;
	}

	public static double T10ShapeFunction(double[][] xyz, double[] L, double[] N, double[][] gradN) {
		double[][] xxyy = MVO.subMatrix(xyz, 0, 3, 0, 2);
		double[][] grad = new double[4][3];
		double volume = gradient3DT4(xxyy, grad);
		volume *= 6.0;
		double[] dxL = MVO.col(grad, 0);
		double[] dyL = MVO.col(grad, 1);
		double[] dzL = MVO.col(grad, 2);

		double L1 = L[0], L2 = L[1], L3 = L[2], L4 = L[3];
		N[0] = L1 * (2 * L1 - 1);
		N[1] = L2 * (2 * L2 - 1);
		N[2] = L3 * (2 * L3 - 1);
		N[3] = L4 * (2 * L4 - 1);
		N[4] = 4 * L1 * L2;
		N[5] = 4 * L2 * L3;
		N[6] = 4 * L3 * L1;
		N[7] = 4 * L1 * L4;
		N[8] = 4 * L2 * L4;
		N[9] = 4 * L3 * L4;
		gradN[0][0] = 4 * dxL[0] * L1 - dxL[0];
		gradN[1][0] = 4 * dxL[1] * L2 - dxL[1];
		gradN[2][0] = 4 * dxL[2] * L3 - dxL[2];
		gradN[3][0] = 4 * dxL[3] * L4 - dxL[3];
		gradN[4][0] = 4 * (dxL[0] * L2 + L1 * dxL[1]);
		gradN[5][0] = 4 * (dxL[1] * L3 + L2 * dxL[2]);
		gradN[6][0] = 4 * (dxL[2] * L1 + L3 * dxL[0]);
		gradN[7][0] = 4 * (dxL[0] * L4 + L1 * dxL[3]);
		gradN[8][0] = 4 * (dxL[1] * L4 + L2 * dxL[3]);
		gradN[9][0] = 4 * (dxL[2] * L4 + L3 * dxL[3]);
		
		gradN[0][1] = 4 * dyL[0] * L1 - dyL[0];
		gradN[1][1] = 4 * dyL[1] * L2 - dyL[1];
		gradN[2][1] = 4 * dyL[2] * L3 - dyL[2];
		gradN[3][1] = 4 * dyL[3] * L4 - dyL[3];
		gradN[4][1] = 4 * (dyL[0] * L2 + L1 * dyL[1]);
		gradN[5][1] = 4 * (dyL[1] * L3 + L2 * dyL[2]);
		gradN[6][1] = 4 * (dyL[2] * L1 + L3 * dyL[0]);
		gradN[7][1] = 4 * (dyL[0] * L4 + L1 * dyL[3]);
		gradN[8][1] = 4 * (dyL[1] * L4 + L2 * dyL[3]);
		gradN[9][1] = 4 * (dyL[2] * L4 + L3 * dyL[3]);
		
		gradN[0][2] = 4 * dzL[0] * L1 - dzL[0];
		gradN[1][2] = 4 * dzL[1] * L2 - dzL[1];
		gradN[2][2] = 4 * dzL[2] * L3 - dzL[2];
		gradN[3][2] = 4 * dzL[3] * L4 - dzL[3];
		gradN[4][2] = 4 * (dzL[0] * L2 + L1 * dzL[1]);
		gradN[5][2] = 4 * (dzL[1] * L3 + L2 * dzL[2]);
		gradN[6][2] = 4 * (dzL[2] * L1 + L3 * dzL[0]);
		gradN[7][2] = 4 * (dzL[0] * L4 + L1 * dzL[3]);
		gradN[8][2] = 4 * (dzL[1] * L4 + L2 * dzL[3]);
		gradN[9][2] = 4 * (dzL[2] * L4 + L3 * dzL[3]);
		return volume;
	}

	public static void L2ShapeFunction(double xi, double[] h, double[] dXi) {
		h[0] = 0.5 * (1 - xi);
		h[1] = 0.5 * (1 + xi);

		dXi[0] = -0.5;
		dXi[1] = 0.5;
	}

	/**
	 * �����ά�ο��ı��εĽ���κ����Լ�����ֵ
	 * 
	 * @param xi
	 * @param et
	 * @param h
	 * @param dXi
	 * @param dEt
	 */
	public static void Q4ShapeFunction(double xi, double et, double[] h, double[] dXi, double[] dEt) {
		h[0] = 0.25 * (1 - xi) * (1 - et);
		h[1] = 0.25 * (1 + xi) * (1 - et);
		h[2] = 0.25 * (1 + xi) * (1 + et);
		h[3] = 0.25 * (1 - xi) * (1 + et);

		dXi[0] = 0.25 * (et - 1);
		dXi[1] = 0.25 * (1 - et);
		dXi[2] = 0.25 * (1 + et);
		dXi[3] = 0.25 * (-1 - et);

		dEt[0] = 0.25 * (xi - 1);
		dEt[1] = 0.25 * (-1 - xi);
		dEt[2] = 0.25 * (1 + xi);
		dEt[3] = 0.25 * (1 - xi);
	}
	
	public static void Q4ShapeFuncionCoef(double[][] xy, double[][] coef) {
//		SRMatrix A = new SRMatrix(4);
//		double[] x = MVO.col(xy, 0);
//		double[] y = MVO.col(xy, 1);
//		for (int i = 0; i < 4; i++) {
//			A.setElementFast(i, 0, 1.0);
//			A.setElementFast(i, 1, x[i]);
//			A.setElementFast(i, 2, y[i]);
//			A.setElementFast(i, 3, x[i] * y[i]);
//		}
//		NewIterSSolver solver = new NewIterSSolver(A);
//		double[] b = new double[4];
//		b[0] = 1.0;
//		solver.CGNE(b, coef[0], 0);
//		b[0] = 0.0;
//		b[1] = 1.0;
//		solver.CGNE(b, coef[1], 0);
//		b[1] = 0.0;
//		b[2] = 1.0;
//		solver.CGNE(b, coef[2], 0);
//		b[2] = 0.0;
//		b[3] = 1.0;
//		solver.CGNE(b, coef[3], 0);
		
		double[][] A = new double[4][];
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		for(int i = 0; i < 4; i++) {
			A[i] = new double[]{1.0, x[i], y[i], x[i] * y[i]};
		}
		
		IterDSolver solver = new IterDSolver();
		double[] b = new double[4];
		b[0] = 1.0;
		solver.CGNE(A, b, coef[0], 0);
		b[0] = 0.0;
		b[1] = 1.0;
		solver.CGNE(A, b, coef[1], 0);
		b[1] = 0.0;
		b[2] = 1.0;
		solver.CGNE(A, b, coef[2], 0);
		b[2] = 0.0;
		b[3] = 1.0;
		solver.CGNE(A, b, coef[3], 0);
	}

	public static void Q4ShapeFunction(double[][] xy, double[] xy0, double[] N) {
//		SRMatrix A = new SRMatrix(4);
//		double[] x = MVO.col(xy, 0);
//		double[] y = MVO.col(xy, 1);
//		for (int i = 0; i < 4; i++) {
//			A.setElementFast(i, 0, 1.0);
//			A.setElementFast(i, 1, x[i]);
//			A.setElementFast(i, 2, y[i]);
//			A.setElementFast(i, 3, x[i] * y[i]);
//		}
//		NewIterSSolver solver = new NewIterSSolver(A);
//		double[] b = new double[4];
//		double[] localXY = { 1.0, xy0[0], xy0[1], xy0[0] * xy0[1] };
//		double[] Nx = new double[4];
//		b[0] = 1.0;
//
//		solver.CGNE(b, Nx, 0);
//		N[0] = MVO.dot_product(localXY, Nx);
//		b[0] = 0.0;
//		b[1] = 1.0;
//		solver.CGNE(b, Nx, 0);
//		N[1] = MVO.dot_product(localXY, Nx);
//		b[1] = 0.0;
//		b[2] = 1.0;
//		solver.CGNE(b, Nx, 0);
//		N[2] = MVO.dot_product(localXY, Nx);
//		b[2] = 0.0;
//		b[3] = 1.0;
//		solver.CGNE(b, Nx, 0);
//		N[3] = MVO.dot_product(localXY, Nx);
		
		double[][] A = new double[4][];
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		for(int i = 0; i < 4; i++) {
			A[i] = new double[]{1.0, x[i], y[i], x[i] * y[i]};
		}
		
		IterDSolver solver = new IterDSolver();
		double[] b = new double[4];
		double[] localXY = { 1.0, xy0[0], xy0[1], xy0[0] * xy0[1] };
		double[] Nx = new double[4];
		b[0] = 1.0;

		solver.CGNE(A, b, Nx, 0);
		N[0] = MVO.dot_product(localXY, Nx);
		b[0] = 0.0;
		b[1] = 1.0;
		solver.CGNE(A, b, Nx, 0);
		N[1] = MVO.dot_product(localXY, Nx);
		b[1] = 0.0;
		b[2] = 1.0;
		solver.CGNE(A, b, Nx, 0);
		N[2] = MVO.dot_product(localXY, Nx);
		b[2] = 0.0;
		b[3] = 1.0;
		solver.CGNE(A, b, Nx, 0);
		N[3] = MVO.dot_product(localXY, Nx);
		
	}

	public static void gradientValue2DQ4(double[][] coord, double[] u0, double[][] gradU) {
		double[] x = MVO.col(coord, 0);
		double[] y = MVO.col(coord, 1);
		double[][] coef = new double[4][4];
		Q4ShapeFuncionCoef(coord, coef);

		double[] dxU = new double[] { coef[0][1] + coef[0][3] * y[0], coef[1][1] + coef[1][3] * y[0],
				coef[2][1] + coef[2][3] * y[0], coef[3][1] + coef[3][3] * y[0] };
		double[] dyU = new double[] { coef[0][2] + coef[0][3] * x[0], coef[1][2] + coef[1][3] * x[0],
				coef[2][2] + coef[2][3] * x[0], coef[3][2] + coef[3][3] * x[0] };
		// (x0, y0)��
		gradU[0][0] = MVO.dot_product(dxU, u0);
		gradU[0][1] = MVO.dot_product(dyU, u0);

		dxU = new double[] { coef[0][1] + coef[0][3] * y[1], coef[1][1] + coef[1][3] * y[1],
				coef[2][1] + coef[2][3] * y[1], coef[3][1] + coef[3][3] * y[1] };
		dyU = new double[] { coef[0][2] + coef[0][3] * x[1], coef[1][2] + coef[1][3] * x[1],
				coef[2][2] + coef[2][3] * x[1], coef[3][2] + coef[3][3] * x[1] };
		gradU[1][0] = MVO.dot_product(dxU, u0);
		gradU[1][1] = MVO.dot_product(dyU, u0);

		dxU = new double[] { coef[0][1] + coef[0][3] * y[2], coef[1][1] + coef[1][3] * y[2],
				coef[2][1] + coef[2][3] * y[2], coef[3][1] + coef[3][3] * y[2] };
		dyU = new double[] { coef[0][2] + coef[0][3] * x[2], coef[1][2] + coef[1][3] * x[2],
				coef[2][2] + coef[2][3] * x[2], coef[3][2] + coef[3][3] * x[2] };
		gradU[2][0] = MVO.dot_product(dxU, u0);
		gradU[2][1] = MVO.dot_product(dyU, u0);

		dxU = new double[] { coef[0][1] + coef[0][3] * y[3], coef[1][1] + coef[1][3] * y[3],
				coef[2][1] + coef[2][3] * y[3], coef[3][1] + coef[3][3] * y[3] };
		dyU = new double[] { coef[0][2] + coef[0][3] * x[3], coef[1][2] + coef[1][3] * x[3],
				coef[2][2] + coef[2][3] * x[3], coef[3][2] + coef[3][3] * x[3] };
		gradU[3][0] = MVO.dot_product(dxU, u0);
		gradU[3][1] = MVO.dot_product(dyU, u0);
	}

	public static void Q9ShapeFunction(double xi, double et, double[] h, double[] dXi, double[] dEt) {
		double center = (1 - xi * xi) * (1 - et * et);
		h[0] = 0.25 * ((1 - xi) * (1 - et) * (-1 - xi - et) + center);
		h[1] = 0.25 * ((1 + xi) * (1 - et) * (-1 + xi - et) + center);
		h[2] = 0.25 * ((1 + xi) * (1 + et) * (-1 + xi + et) + center);
		h[3] = 0.25 * ((1 - xi) * (1 + et) * (-1 - xi + et) + center);
		h[4] = 0.5 * ((1 - xi * xi) * (1 - et) - center);
		h[5] = 0.5 * ((1 + xi) * (1 - et * et) - center);
		h[6] = 0.5 * ((1 - xi * xi) * (1 + et) - center);
		h[7] = 0.5 * ((1 - xi) * (1 - et * et) - center);
		h[8] = center;

		dXi[0] = 0.25 * (1 - et) * (2 * xi + et) - 0.5 * xi * (1 - et * et);
		dXi[1] = 0.25 * (1 - et) * (2 * xi - et) - 0.5 * xi * (1 - et * et);
		dXi[2] = 0.25 * (1 + et) * (2 * xi + et) - 0.5 * xi * (1 - et * et);
		dXi[3] = 0.25 * (1 + et) * (2 * xi - et) - 0.5 * xi * (1 - et * et);

		dXi[4] = xi * (et - 1)        + xi * (1 - et * et);
		dXi[5] = 0.5 * (1 - et * et)  + xi * (1 - et * et);
		dXi[6] = xi * (-et - 1)       + xi * (1 - et * et);
		dXi[7] = 0.5 * (-1 + et * et) + xi * (1 - et * et);
		dXi[8] = -2 * xi * (1 - et * et);

		dEt[0] = 0.25 * (1 - xi) * (xi + 2 * et)  - 0.5 * et * (1 - xi * xi);
		dEt[1] = 0.25 * (1 + xi) * (-xi + 2 * et) - 0.5 * et * (1 - xi * xi);
		dEt[2] = 0.25 * (1 + xi) * (xi + 2 * et)  - 0.5 * et * (1 - xi * xi);
		dEt[3] = 0.25 * (1 - xi) * (-xi + 2 * et) - 0.5 * et * (1 - xi * xi);

		dEt[4] = 0.5 * (-1 + xi * xi) + et * (1 - xi * xi);
		dEt[5] = et * (-xi - 1)       + et * (1 - xi * xi);
		dEt[6] = 0.5 * (1 - xi * xi)  + et * (1 - xi * xi);
		dEt[7] = et * (xi - 1)        + et * (1 - xi * xi);
		dEt[8] = -2 * et * (1 - xi * xi);
	}

	public static void Q8ShapeFunction(double xi, double et, double[] h, double[] dXi, double[] dEt) {
		h[0] = 0.25 * (1 - xi) * (1 - et) * (-1 - xi - et);
		h[1] = 0.25 * (1 + xi) * (1 - et) * (-1 + xi - et);
		h[2] = 0.25 * (1 + xi) * (1 + et) * (-1 + xi + et);
		h[3] = 0.25 * (1 - xi) * (1 + et) * (-1 - xi + et);
		h[4] = 0.5 * (1 - xi * xi) * (1 - et);
		h[5] = 0.5 * (1 + xi) * (1 - et * et);
		h[6] = 0.5 * (1 - xi * xi) * (1 + et);
		h[7] = 0.5 * (1 - xi) * (1 - et * et);

		dXi[0] = 0.25 * (1 - et) * (2 * xi + et);
		dXi[1] = 0.25 * (1 - et) * (2 * xi - et);
		dXi[2] = 0.25 * (1 + et) * (2 * xi + et);
		dXi[3] = 0.25 * (1 + et) * (2 * xi - et);

		dXi[4] = xi * (et - 1);
		dXi[5] = 0.5 * (1 - et * et);
		dXi[6] = xi * (-et - 1);
		dXi[7] = 0.5 * (-1 + et * et);

		dEt[0] = 0.25 * (1 - xi) * (xi + 2 * et);
		dEt[1] = 0.25 * (1 + xi) * (-xi + 2 * et);
		dEt[2] = 0.25 * (1 + xi) * (xi + 2 * et);
		dEt[3] = 0.25 * (1 - xi) * (-xi + 2 * et);

		dEt[4] = 0.5 * (-1 + xi * xi);
		dEt[5] = et * (-xi - 1);
		dEt[6] = 0.5 * (1 - xi * xi);
		dEt[7] = et * (xi - 1);
	}

	public static void Q8ShapeFunctionCoef(double[][] xy, double[][] coef) {
//		SRMatrix A = new SRMatrix(8);
//		double[] x = MVO.col(xy, 0);
//		double[] y = MVO.col(xy, 1);
//		for (int i = 0; i < 8; i++) {
//			A.setElementFast(i, 0, 1.0);
//			A.setElementFast(i, 1, x[i]);
//			A.setElementFast(i, 2, y[i]);
//			A.setElementFast(i, 3, x[i] * x[i]);
//			A.setElementFast(i, 4, x[i] * y[i]);
//			A.setElementFast(i, 5, y[i] * y[i]);
//			A.setElementFast(i, 6, x[i] * x[i] * y[i]);
//			A.setElementFast(i, 7, x[i] * y[i] * y[i]);
//		}
//
//		NewIterSSolver solver = new NewIterSSolver(A);
//		double[] b = new double[8];
//		b[0] = 1.0;
//		for (int i = 0; i < 7; i++) {
//			solver.CGNE(b, coef[i], 0);
//			b[i] = 0.0;
//			b[i + 1] = 1.0;
//		}
//		solver.CGNE(b, coef[7], 0);
		
		double[][] A = new double[8][];
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		for(int i = 0; i < 8; i++) {
			A[i] = new double[] {1.0, x[i], y[i], x[i] * x[i], x[i] * y[i], y[i] * y[i], x[i] * x[i] * y[i], x[i] * y[i] * y[i]};
		}
		IterDSolver solver = new IterDSolver();
//		double[] b = new double[8];
//		b[0] = 1.0;
//		for (int i = 0; i < 7; i++) {
//			solver.CGNE(A, b, coef[i], 0);
//			b[i] = 0.0;
//			b[i + 1] = 1.0;
//		}
//		solver.CGNE(A, b, coef[7], 0);		
		
		double[] b = new double[8];
		b[0] = 1.0;
		for (int i = 0; i < 7; i++) {
			solver.CGNE(A, b, coef[i], 0);
			b[i] = 0.0;
			b[i + 1] = 1.0;
		}
		solver.CGNE(A, b, coef[7], 0);
	}

	public static void Q8ShapeFunction(double[][] xy, double[] xy0, double[] N) {
		double[][] coef = new double[8][8];
		Q8ShapeFunctionCoef(xy, coef);
		double x = xy0[0], y = xy0[1];
		double[] v = new double[] { 1.0, x, y, x * x, x * y, y * y, x * x * y, x * y * y };
		MVO.matmul(coef, v, N);
	}

	public static void gradientValue2DQ8(double[][] coord, double[] u0, double[][] gradU) {
		double[][] coef = new double[8][8];
		Q8ShapeFunctionCoef(coord, coef);
		double x, y;
		double[] dxU = new double[8];
		double[] dyU = new double[8];
		for (int i = 0; i < 8; i++) {
			x = coord[i][0];
			y = coord[i][1];
			for (int j = 0; j < 8; j++) {
				dxU[j] = coef[j][1] + 2 * coef[j][3] * x + coef[j][4] * y + 2 * coef[j][6] * x * y + coef[j][7] * y * y;
				dyU[j] = coef[j][2] + coef[j][4] * x + 2 * coef[j][5] * y + coef[j][6] * x * x + 2 * coef[j][7] * x * y;
			}
			gradU[i][0] = MVO.dot_product(dxU, u0);
			gradU[i][1] = MVO.dot_product(dyU, u0);
		}
	}

	public static void secondDerivativeValue2DQ8(double[][] coord, double[] u0, double[][] sDU) {
		double[][] coef = new double[8][8];
		Q8ShapeFunctionCoef(coord, coef);
		double x, y;
		double[] dxxU = new double[8];
		double[] dxyU = new double[8];
		double[] dyyU = new double[8];
		for (int i = 0; i < 8; i++) {
			x = coord[i][0];
			y = coord[i][1];
			for (int j = 0; j < 8; j++) {
				dxxU[j] = 2 * coef[j][3] + 2 * coef[j][6] * y;
				dxyU[j] = coef[j][4] + 2 * coef[j][6] * x + 2 * coef[j][7] * y;
				dyyU[j] = 2 * coef[j][5] + 2 * coef[j][7] * x;

			}
			sDU[i][0] = MVO.dot_product(dxxU, u0);
			sDU[i][1] = MVO.dot_product(dxyU, u0);
			sDU[i][2] = MVO.dot_product(dyyU, u0);
		}
	}
	
	public static void firstAndSecondDerivativeValue2DQ8(double[][] coord, double[] u0, double[][] fDU, double[][] sDU) {
		double[][] coef = new double[8][8];
		Q8ShapeFunctionCoef(coord, coef);
		double x, y;
		double[] dxU = new double[8];
		double[] dyU = new double[8];
		double[] dxxU = new double[8];
		double[] dxyU = new double[8];
		double[] dyyU = new double[8];
		for (int i = 0; i < 8; i++) {
			x = coord[i][0];
			y = coord[i][1];
			for (int j = 0; j < 8; j++) {
				dxU[j] = coef[j][1] + 2 * coef[j][3] * x + coef[j][4] * y + 2 * coef[j][6] * x * y + coef[j][7] * y * y;
				dyU[j] = coef[j][2] + coef[j][4] * x + 2 * coef[j][5] * y + coef[j][6] * x * x + 2 * coef[j][7] * x * y;
				dxxU[j] = 2 * coef[j][3] + 2 * coef[j][6] * y;
				dxyU[j] = coef[j][4] + 2 * coef[j][6] * x + 2 * coef[j][7] * y;
				dyyU[j] = 2 * coef[j][5] + 2 * coef[j][7] * x;
			}
			fDU[i][0] = MVO.dot_product(dxU, u0);
			fDU[i][1] = MVO.dot_product(dyU, u0);
			sDU[i][0] = MVO.dot_product(dxxU, u0);
			sDU[i][1] = MVO.dot_product(dxyU, u0);
			sDU[i][2] = MVO.dot_product(dyyU, u0);
		}
	}

	/**
	 * 
	 * @param gradXY
	 *            6����������������κ����ݶ�ֵ
	 * @param u0
	 * @param gradU
	 */
	public static void gradientValue2DT6(double[][] gradXY, double[] u0, double[][] gradU) {
		double[] dx = MVO.col(gradXY, 0);
		double[] dy = MVO.col(gradXY, 1);

		double[] dxN = new double[] { 3 * dx[0], -dx[1], -dx[2], 4 * dx[1], 0, 4 * dx[2] };
		double[] dyN = new double[] { 3 * dy[0], -dy[1], -dy[2], 4 * dy[1], 0, 4 * dy[2] };

		gradU[0][0] = MVO.dot_product(u0, dxN);
		gradU[0][1] = MVO.dot_product(u0, dyN);

		dxN = new double[] { -dx[0], 3 * dx[1], -dx[2], 4 * dx[0], 4 * dx[2], 0 };
		dyN = new double[] { -dy[0], 3 * dy[1], -dy[2], 4 * dy[0], 4 * dy[2], 0 };

		gradU[1][0] = MVO.dot_product(u0, dxN);
		gradU[1][1] = MVO.dot_product(u0, dyN);

		dxN = new double[] { -dx[0], -dx[1], 3 * dx[2], 0, 4 * dx[1], 4 * dx[0] };
		dyN = new double[] { -dy[0], -dy[1], 3 * dy[2], 0, 4 * dy[1], 4 * dy[0] };

		gradU[2][0] = MVO.dot_product(u0, dxN);
		gradU[2][1] = MVO.dot_product(u0, dyN);

		dxN = new double[] { dx[0], dx[1], -dx[2], 2 * (dx[0] + dx[1]), 2 * dx[2], 2 * dx[2] };
		dyN = new double[] { dy[0], dy[1], -dy[2], 2 * (dy[0] + dy[1]), 2 * dy[2], 2 * dy[2] };
		gradU[3][0] = MVO.dot_product(u0, dxN);
		gradU[3][1] = MVO.dot_product(u0, dyN);

		dxN = new double[] { -dx[0], dx[1], dx[2], 2 * dx[0], 2 * (dx[1] + dx[2]), 2 * dx[0] };
		dyN = new double[] { -dy[0], dy[1], dy[2], 2 * dy[0], 2 * (dy[1] + dy[2]), 2 * dy[0] };
		gradU[4][0] = MVO.dot_product(u0, dxN);
		gradU[4][1] = MVO.dot_product(u0, dyN);

		dxN = new double[] { dx[0], -dx[1], dx[2], 2 * dx[1], 2 * dx[1], 2 * (dx[0] + dx[2]) };
		dyN = new double[] { dy[0], -dy[1], dy[2], 2 * dy[1], 2 * dy[1], 2 * (dy[0] + dy[2]) };
		gradU[5][0] = MVO.dot_product(u0, dxN);
		gradU[5][1] = MVO.dot_product(u0, dyN);
	}

	public static void secondDerivative2DT6(double[][] gradXY, double[] u0, double[] sdU) {
		double[] dx = MVO.col(gradXY, 0);
		double[] dy = MVO.col(gradXY, 1);

		double[] dxx = { 4 * dx[0] * dx[0], 4 * dx[1] * dx[1], 4 * dx[2] * dx[2], 8 * dx[0] * dx[1], 8 * dx[1] * dx[2],
				8 * dx[0] * dx[2] };
		double[] dyy = { 4 * dy[0] * dy[0], 4 * dy[1] * dy[1], 4 * dy[2] * dy[2], 8 * dy[0] * dy[1], 8 * dy[1] * dy[2],
				8 * dy[0] * dy[2] };
		double[] dxy = { 4 * dx[0] * dy[0], 4 * dx[1] * dy[1], 4 * dx[2] * dy[2], 4 * (dx[0] * dy[1] + dy[0] * dx[1]),
				4 * (dx[1] * dy[2] + dy[1] * dx[2]), 4 * (dx[0] * dy[2] + dy[0] * dx[2]) };
		sdU[0] = MVO.dot_product(u0, dxx); // 11
		sdU[1] = MVO.dot_product(u0, dyy); // 22
		sdU[2] = MVO.dot_product(u0, dxy); // 12
	}
	
	public static void gradientValue3DT10(double[][] gradXYZ, double[] u0, double[][] gradU) {
		double[] dx = MVO.col(gradXYZ, 0);
		double[] dy = MVO.col(gradXYZ, 1);
		double[] dz = MVO.col(gradXYZ, 2);
		
		//N1 = 1, N2 = 0, N3 = 0, N4 = 0
		double[] dxN = new double[] { 3 * dx[0], -dx[1], -dx[2], -dx[3], 4 * dx[1], 0, 4 * dx[2], 4 * dx[3], 0, 0};
		double[] dyN = new double[] { 3 * dy[0], -dy[1], -dy[2], -dy[3], 4 * dy[1], 0, 4 * dy[2], 4 * dy[3], 0, 0};
		double[] dzN = new double[] { 3 * dz[0], -dz[1], -dz[2], -dz[3], 4 * dz[1], 0, 4 * dz[2], 4 * dz[3], 0, 0};

		gradU[0][0] = MVO.dot_product(u0, dxN);
		gradU[0][1] = MVO.dot_product(u0, dyN);
		gradU[0][2] = MVO.dot_product(u0, dzN);
		
		//N1 = 0, N2 = 1, N3 = 0, N4 = 0
		dxN = new double[] {-dx[0], 3 * dx[1], -dx[2], -dx[3], 4 * dx[0], 4 * dx[2], 0, 0, 4 * dx[3], 0};
		dyN = new double[] {-dy[0], 3 * dy[1], -dy[2], -dy[3], 4 * dy[0], 4 * dy[2], 0, 0, 4 * dy[3], 0};
		dzN = new double[] {-dz[0], 3 * dz[1], -dz[2], -dz[3], 4 * dz[0], 4 * dz[2], 0, 0, 4 * dz[3], 0};
	
		gradU[1][0] = MVO.dot_product(u0, dxN);
		gradU[1][1] = MVO.dot_product(u0, dyN);
		gradU[1][2] = MVO.dot_product(u0, dzN);
		
		//N1 = 0, N2 = 0, N3 = 1, N4 = 0
		dxN = new double[] {-dx[0], -dx[1], 3 * dx[2], -dx[3], 0, 4 * dx[1], 4 * dx[0], 0, 0, 4 * dx[3]};
		dyN = new double[] {-dy[0], -dy[1], 3 * dy[2], -dy[3], 0, 4 * dy[1], 4 * dy[0], 0, 0, 4 * dy[3]};
		dzN = new double[] {-dz[0], -dz[1], 3 * dz[2], -dz[3], 0, 4 * dz[1], 4 * dz[0], 0, 0, 4 * dz[3]};
		
		gradU[2][0] = MVO.dot_product(u0, dxN);
		gradU[2][1] = MVO.dot_product(u0, dyN);
		gradU[2][2] = MVO.dot_product(u0, dzN);
		
		//N1 = 0, N2 = 0, N3 = 0, N4 = 1
		dxN = new double[] {-dx[0], -dx[1], -dx[2], 3 * dx[3], 0, 0, 0, 4 * dx[0], 4 * dx[1], 4 * dx[2]};
		dyN = new double[] {-dy[0], -dy[1], -dy[2], 3 * dy[3], 0, 0, 0, 4 * dy[0], 4 * dy[1], 4 * dy[2]};
		dzN = new double[] {-dz[0], -dz[1], -dz[2], 3 * dz[3], 0, 0, 0, 4 * dz[0], 4 * dz[1], 4 * dz[2]};
		
		gradU[3][0] = MVO.dot_product(u0, dxN);
		gradU[3][1] = MVO.dot_product(u0, dyN);
		gradU[3][2] = MVO.dot_product(u0, dzN);
		
		//N1 = 0.5, N2 = 0.5, N3 = 0, N4 = 0
		dxN = new double[] {dx[0], dx[1], -dx[2], -dx[3], 2 * (dx[0] + dx[1]), 2 * dx[2], 2 * dx[2], 2 * dx[3], 2 * dx[3], 0};
		dyN = new double[] {dy[0], dy[1], -dy[2], -dy[3], 2 * (dy[0] + dy[1]), 2 * dy[2], 2 * dy[2], 2 * dy[3], 2 * dy[3], 0};
		dzN = new double[] {dz[0], dz[1], -dz[2], -dz[3], 2 * (dz[0] + dz[1]), 2 * dz[2], 2 * dz[2], 2 * dz[3], 2 * dz[3], 0};
		
		gradU[4][0] = MVO.dot_product(u0, dxN);
		gradU[4][1] = MVO.dot_product(u0, dyN);
		gradU[4][2] = MVO.dot_product(u0, dzN);
		
		//N1 = 0, N2 = 0.5, N3 = 0.5, N4 = 0
		dxN = new double[] {-dx[0], dx[1], dx[2], -dx[3], 2 * dx[0], 2 * (dx[1] + dx[2]), 2 * dx[0], 0, 2 * dx[3], 2 * dx[3]};
		dyN = new double[] {-dy[0], dy[1], dy[2], -dy[3], 2 * dy[0], 2 * (dy[1] + dy[2]), 2 * dy[0], 0, 2 * dy[3], 2 * dy[3]};
		dzN = new double[] {-dz[0], dz[1], dz[2], -dz[3], 2 * dz[0], 2 * (dz[1] + dz[2]), 2 * dz[0], 0, 2 * dz[3], 2 * dz[3]};

		gradU[5][0] = MVO.dot_product(u0, dxN);
		gradU[5][1] = MVO.dot_product(u0, dyN);
		gradU[5][2] = MVO.dot_product(u0, dzN);
		
		//N1 = 0.5, N2 = 0, N3 = 0.5, N4 = 0
		dxN = new double[] {dx[0], -dx[1], dx[2], -dx[3], 2 * dx[1], 2 * dx[1], 2 * (dx[0] + dx[2]), 2 * dx[3], 0, 2 * dx[3]};
		dyN = new double[] {dy[0], -dy[1], dy[2], -dy[3], 2 * dy[1], 2 * dy[1], 2 * (dy[0] + dy[2]), 2 * dy[3], 0, 2 * dy[3]};
		dzN = new double[] {dz[0], -dz[1], dz[2], -dz[3], 2 * dz[1], 2 * dz[1], 2 * (dz[0] + dz[2]), 2 * dz[3], 0, 2 * dz[3]};
		
		gradU[6][0] = MVO.dot_product(u0, dxN);
		gradU[6][1] = MVO.dot_product(u0, dyN);
		gradU[6][2] = MVO.dot_product(u0, dzN);		
		
		//N1 = 0.5, N2 = 0, N3 = 0, N4 = 0.5
		dxN = new double[] {dx[0], -dx[1], -dx[2], dx[3], 2 * dx[1], 0, 2 * dx[2], 2 * (dx[0] + dx[3]), 2 * dx[1], 2 * dx[2]};
		dyN = new double[] {dy[0], -dy[1], -dy[2], dy[3], 2 * dy[1], 0, 2 * dy[2], 2 * (dy[0] + dy[3]), 2 * dy[1], 2 * dy[2]};
		dzN = new double[] {dz[0], -dz[1], -dz[2], dz[3], 2 * dz[1], 0, 2 * dz[2], 2 * (dz[0] + dz[3]), 2 * dz[1], 2 * dz[2]};

		gradU[7][0] = MVO.dot_product(u0, dxN);
		gradU[7][1] = MVO.dot_product(u0, dyN);
		gradU[7][2] = MVO.dot_product(u0, dzN);
		
		//N1 = 0, N2 = 0.5, N3 = 0, N4 = 0.5
		dxN = new double[] {-dx[0], dx[1], -dx[2], dx[3], 2 * dx[0], 2 * dx[2], 0, 2 * dx[0], 2 * (dx[1] + dx[3]), 2 * dx[2]};
		dyN = new double[] {-dy[0], dy[1], -dy[2], dy[3], 2 * dy[0], 2 * dy[2], 0, 2 * dy[0], 2 * (dy[1] + dy[3]), 2 * dy[2]};
		dzN = new double[] {-dz[0], dz[1], -dz[2], dz[3], 2 * dz[0], 2 * dz[2], 0, 2 * dz[0], 2 * (dz[1] + dz[3]), 2 * dz[2]};
		
		gradU[8][0] = MVO.dot_product(u0, dxN);
		gradU[8][1] = MVO.dot_product(u0, dyN);
		gradU[8][2] = MVO.dot_product(u0, dzN);
		
		//N1 = 0, N2 = 0, N3 = 0.5, N4 = 0.5
		dxN = new double[] {-dx[0], -dx[1], dx[2], dx[3], 0, 2 * dx[1], 2 * dx[0], 2 * dx[0], 2 * dx[1], 2 * (dx[2] + dx[3])};
		dyN = new double[] {-dy[0], -dy[1], dy[2], dy[3], 0, 2 * dy[1], 2 * dy[0], 2 * dy[0], 2 * dy[1], 2 * (dy[2] + dy[3])};
		dzN = new double[] {-dz[0], -dz[1], dz[2], dz[3], 0, 2 * dz[1], 2 * dz[0], 2 * dz[0], 2 * dz[1], 2 * (dz[2] + dz[3])};
		
		gradU[9][0] = MVO.dot_product(u0, dxN);
		gradU[9][1] = MVO.dot_product(u0, dyN);
		gradU[9][2] = MVO.dot_product(u0, dzN);
		
	}
	
	public static void secondDerivative3DT10(double[][] gradXYZ, double[] u0, double[] sdU) {
		double[] dx = MVO.col(gradXYZ, 0);
		double[] dy = MVO.col(gradXYZ, 1);
		double[] dz = MVO.col(gradXYZ, 2);

		double[] dxx = { 4 * dx[0] * dx[0], 4 * dx[1] * dx[1], 4 * dx[2] * dx[2], 4 * dx[3] * dx[3], 
				8 * dx[0] * dx[1], 8 * dx[1] * dx[2], 8 * dx[0] * dx[2], 8 * dx[0] * dx[3], 8 * dx[1] * dx[3], 8 * dx[2] * dx[3]};
		double[] dyy = { 4 * dy[0] * dy[0], 4 * dy[1] * dy[1], 4 * dy[2] * dy[2], 4 * dy[3] * dy[3], 
				8 * dy[0] * dy[1], 8 * dy[1] * dy[2], 8 * dy[0] * dy[2], 8 * dy[0] * dy[3], 8 * dy[1] * dy[3], 8 * dy[2] * dy[3]};
		double[] dzz = { 4 * dz[0] * dz[0], 4 * dz[1] * dz[1], 4 * dz[2] * dz[2], 4 * dz[3] * dz[3], 
				8 * dz[0] * dz[1], 8 * dz[1] * dz[2], 8 * dz[0] * dz[2], 8 * dz[0] * dz[3], 8 * dz[1] * dz[3], 8 * dz[2] * dz[3]};
		
		double[] dxy = { 4 * dx[0] * dy[0], 4 * dx[1] * dy[1], 4 * dx[2] * dy[2], 4 * dx[3] * dy[3], 
				4 * (dx[0] * dy[1] + dy[0] * dx[1]), 
				4 * (dx[1] * dy[2] + dy[1] * dx[2]), 
				4 * (dx[0] * dy[2] + dy[0] * dx[2]), 
				4 * (dx[0] * dy[3] + dy[0] * dx[3]), 
				4 * (dx[1] * dy[3] + dy[1] * dx[3]), 
				4 * (dx[2] * dy[3] + dy[2] * dx[3])};
		double[] dyz = { 4 * dz[0] * dy[0], 4 * dz[1] * dy[1], 4 * dz[2] * dy[2], 4 * dz[3] * dy[3], 
				4 * (dz[0] * dy[1] + dy[0] * dz[1]), 
				4 * (dz[1] * dy[2] + dy[1] * dz[2]), 
				4 * (dz[0] * dy[2] + dy[0] * dz[2]), 
				4 * (dz[0] * dy[3] + dy[0] * dz[3]), 
				4 * (dz[1] * dy[3] + dy[1] * dz[3]), 
				4 * (dz[2] * dy[3] + dy[2] * dz[3])};
		double[] dxz = { 4 * dz[0] * dx[0], 4 * dz[1] * dx[1], 4 * dz[2] * dx[2], 4 * dz[3] * dx[3], 
				4 * (dz[0] * dx[1] + dx[0] * dz[1]), 
				4 * (dz[1] * dx[2] + dx[1] * dz[2]), 
				4 * (dz[0] * dx[2] + dx[0] * dz[2]), 
				4 * (dz[0] * dx[3] + dx[0] * dz[3]), 
				4 * (dz[1] * dx[3] + dx[1] * dz[3]), 
				4 * (dz[2] * dx[3] + dx[2] * dz[3])};
		sdU[0] = MVO.dot_product(u0, dxx); // 11
		sdU[1] = MVO.dot_product(u0, dyy); // 22
		sdU[2] = MVO.dot_product(u0, dzz); // 33
		sdU[3] = MVO.dot_product(u0, dxy); // 12
		sdU[4] = MVO.dot_product(u0, dyz); // 23
		sdU[5] = MVO.dot_product(u0, dxz); // 13
	}

	public static void gradientValue3DH8(double[][] coord, double[] eleU, double[][] grad) {
		double[][] coef = new double[8][8];
		H8ShapeFunctionCoef(coord, coef);
		double x, y, z;
		double[] dxU = new double[8];
		double[] dyU = new double[8];
		double[] dzU = new double[8];
		for (int i = 0; i < 8; i++) {
			x = coord[i][0];
			y = coord[i][1];
			z = coord[i][2];
			for (int j = 0; j < 8; j++) {
				dxU[j] = coef[j][1] + coef[j][4] * y + coef[j][6] * z + coef[j][7] * y * z;
				dyU[j] = coef[j][2] + coef[j][4] * x + coef[j][5] * z + coef[j][7] * x * z;
				dzU[j] = coef[j][3] + coef[j][5] * y + coef[j][6] * x + coef[j][7] * x * y;
			}
			grad[i][0] = MVO.dot_product(dxU, eleU);
			grad[i][1] = MVO.dot_product(dyU, eleU);
			grad[i][2] = MVO.dot_product(dzU, eleU);
		}
	}

	public static void H8ShapeFunction(double xi, double et, double zt, double[] h, double[] dXi, double[] dEt,
			double[] dZt) {
		h[0] = 0.125 * (1 - xi) * (1 - et) * (1 - zt);
		h[1] = 0.125 * (1 + xi) * (1 - et) * (1 - zt);
		h[2] = 0.125 * (1 + xi) * (1 + et) * (1 - zt);
		h[3] = 0.125 * (1 - xi) * (1 + et) * (1 - zt);

		h[4] = 0.125 * (1 - xi) * (1 - et) * (1 + zt);
		h[5] = 0.125 * (1 + xi) * (1 - et) * (1 + zt);
		h[6] = 0.125 * (1 + xi) * (1 + et) * (1 + zt);
		h[7] = 0.125 * (1 - xi) * (1 + et) * (1 + zt);

		dXi[0] = 0.125 * (et - 1) * (1 - zt);
		dXi[1] = 0.125 * (1 - et) * (1 - zt);
		dXi[2] = 0.125 * (1 + et) * (1 - zt);
		dXi[3] = 0.125 * (-1 - et) * (1 - zt);
		dXi[4] = 0.125 * (et - 1) * (1 + zt);
		dXi[5] = 0.125 * (1 - et) * (1 + zt);
		dXi[6] = 0.125 * (1 + et) * (1 + zt);
		dXi[7] = 0.125 * (-1 - et) * (1 + zt);

		dEt[0] = 0.125 * (xi - 1) * (1 - zt);
		dEt[1] = 0.125 * (-1 - xi) * (1 - zt);
		dEt[2] = 0.125 * (1 + xi) * (1 - zt);
		dEt[3] = 0.125 * (1 - xi) * (1 - zt);
		dEt[4] = 0.125 * (xi - 1) * (1 + zt);
		dEt[5] = 0.125 * (-1 - xi) * (1 + zt);
		dEt[6] = 0.125 * (1 + xi) * (1 + zt);
		dEt[7] = 0.125 * (1 - xi) * (1 + zt);

		dZt[0] = -0.125 * (1 - xi) * (1 - et);
		dZt[1] = -0.125 * (1 + xi) * (1 - et);
		dZt[2] = -0.125 * (1 + xi) * (1 + et);
		dZt[3] = -0.125 * (1 - xi) * (1 + et);
		dZt[4] = 0.125 * (1 - xi) * (1 - et);
		dZt[5] = 0.125 * (1 + xi) * (1 - et);
		dZt[6] = 0.125 * (1 + xi) * (1 + et);
		dZt[7] = 0.125 * (1 - xi) * (1 + et);
	}

	public static void H8ShapeFunctionCoef(double[][] xy, double[][] coef) {
		double[][] A = new double[8][];
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] z = MVO.col(xy, 2);
		for (int i = 0; i < 8; i++) {
			A[i] = new double[] {1.0, x[i], y[i], z[i], x[i] * y[i], y[i] * z[i], x[i] * z[i], x[i] * y[i] * z[i]};
		}

		IterDSolver solver = new IterDSolver();
		double[] b = new double[8];
		b[0] = 1.0;
		for (int i = 0; i < 7; i++) {
			solver.CGNE(A, b, coef[i], 0);
			b[i] = 0.0;
			b[i + 1] = 1.0;
		}
		solver.CGNE(A, b, coef[7], 0);
	}

	public static void H8ShapeFunction(double[][] xyz, double[] xyz0, double[] N) {
		double[][] coef = new double[8][8];
		H8ShapeFunctionCoef(xyz, coef);
		double x = xyz0[0], y = xyz0[1], z = xyz0[2];
		double[] v = new double[] { 1.0, x, y, z, x * y, y * z, z * x, x * y * z };
		MVO.matmul(coef, v, N);
	}


	public static void H20ShapeFunction(double xi, double et, double zt, double[] h, double[] dXi, double[] dEt,
			double[] dZt) {
		h[0] = 0.125 * (1 - xi) * (1 - et) * (1 - zt) * (-2 - xi - et - zt);
		h[1] = 0.125 * (1 + xi) * (1 - et) * (1 - zt) * (-2 + xi - et - zt);
		h[2] = 0.125 * (1 + xi) * (1 + et) * (1 - zt) * (-2 + xi + et - zt);
		h[3] = 0.125 * (1 - xi) * (1 + et) * (1 - zt) * (-2 - xi + et - zt);

		h[4] = 0.125 * (1 - xi) * (1 - et) * (1 + zt) * (-2 - xi - et + zt);
		h[5] = 0.125 * (1 + xi) * (1 - et) * (1 + zt) * (-2 + xi - et + zt);
		h[6] = 0.125 * (1 + xi) * (1 + et) * (1 + zt) * (-2 + xi + et + zt);
		h[7] = 0.125 * (1 - xi) * (1 + et) * (1 + zt) * (-2 - xi + et + zt);

		h[8] = 0.25 * (1 - xi * xi) * (1 - et) * (1 - zt);
		h[10] = 0.25 * (1 - xi * xi) * (1 + et) * (1 - zt);
		h[14] = 0.25 * (1 - xi * xi) * (1 + et) * (1 + zt);
		h[12] = 0.25 * (1 - xi * xi) * (1 - et) * (1 + zt);

		h[9] = 0.25 * (1 - et * et) * (1 + xi) * (1 - zt);
		h[13] = 0.25 * (1 - et * et) * (1 + xi) * (1 + zt);
		h[15] = 0.25 * (1 - et * et) * (1 - xi) * (1 + zt);
		h[11] = 0.25 * (1 - et * et) * (1 - xi) * (1 - zt);

		h[16] = 0.25 * (1 - zt * zt) * (1 - xi) * (1 - et);
		h[17] = 0.25 * (1 - zt * zt) * (1 + xi) * (1 - et);
		h[18] = 0.25 * (1 - zt * zt) * (1 + xi) * (1 + et);
		h[19] = 0.25 * (1 - zt * zt) * (1 - xi) * (1 + et);

		dXi[0] = -0.125 * (1 - et) * (1 - zt) * (-1 - 2 * xi - et - zt);
		dXi[1] = 0.125 * (1 - et) * (1 - zt) * (-1 + 2 * xi - et - zt);
		dXi[2] = 0.125 * (1 + et) * (1 - zt) * (-1 + 2 * xi + et - zt);
		dXi[3] = -0.125 * (1 + et) * (1 - zt) * (-1 - 2 * xi + et - zt);

		dXi[4] = -0.125 * (1 - et) * (1 + zt) * (-1 - 2 * xi - et + zt);
		dXi[5] = 0.125 * (1 - et) * (1 + zt) * (-1 + 2 * xi - et + zt);
		dXi[6] = 0.125 * (1 + et) * (1 + zt) * (-1 + 2 * xi + et + zt);
		dXi[7] = -0.125 * (1 + et) * (1 + zt) * (-1 - 2 * xi + et + zt);

		dXi[8] = -0.5 * xi * (1 - et) * (1 - zt);
		dXi[10] = -0.5 * xi * (1 + et) * (1 - zt);
		dXi[14] = -0.5 * xi * (1 + et) * (1 + zt);
		dXi[12] = -0.5 * xi * (1 - et) * (1 + zt);

		dXi[9] = 0.25 * (1 - et * et) * (1 - zt);
		dXi[13] = 0.25 * (1 - et * et) * (1 + zt);
		dXi[15] = -0.25 * (1 - et * et) * (1 + zt);
		dXi[11] = -0.25 * (1 - et * et) * (1 - zt);

		dXi[16] = -0.25 * (1 - zt * zt) * (1 - et);
		dXi[17] = 0.25 * (1 - zt * zt) * (1 - et);
		dXi[18] = 0.25 * (1 - zt * zt) * (1 + et);
		dXi[19] = -0.25 * (1 - zt * zt) * (1 + et);

		dEt[0] = -0.125 * (1 - xi) * (1 - zt) * (-1 - xi - 2 * et - zt);
		dEt[1] = -0.125 * (1 + xi) * (1 - zt) * (-1 + xi - 2 * et - zt);
		dEt[2] = 0.125 * (1 + xi) * (1 - zt) * (-1 + xi + 2 * et - zt);
		dEt[3] = 0.125 * (1 - xi) * (1 - zt) * (-1 - xi + 2 * et - zt);

		dEt[4] = -0.125 * (1 - xi) * (1 + zt) * (-1 - xi - 2 * et + zt);
		dEt[5] = -0.125 * (1 + xi) * (1 + zt) * (-1 + xi - 2 * et + zt);
		dEt[6] = 0.125 * (1 + xi) * (1 + zt) * (-1 + xi + 2 * et + zt);
		dEt[7] = 0.125 * (1 - xi) * (1 + zt) * (-1 - xi + 2 * et + zt);

		dEt[8] = -0.25 * (1 - xi * xi) * (1 - zt);
		dEt[10] = 0.25 * (1 - xi * xi) * (1 - zt);
		dEt[14] = 0.25 * (1 - xi * xi) * (1 + zt);
		dEt[12] = -0.25 * (1 - xi * xi) * (1 + zt);

		dEt[9] = -0.5 * et * (1 + xi) * (1 - zt);
		dEt[13] = -0.5 * et * (1 + xi) * (1 + zt);
		dEt[15] = -0.5 * et * (1 - xi) * (1 + zt);
		dEt[11] = -0.5 * et * (1 - xi) * (1 - zt);

		dEt[16] = -0.25 * (1 - zt * zt) * (1 - xi);
		dEt[17] = -0.25 * (1 - zt * zt) * (1 + xi);
		dEt[18] = 0.25 * (1 - zt * zt) * (1 + xi);
		dEt[19] = 0.25 * (1 - zt * zt) * (1 - xi);

		dZt[0] = -0.125 * (1 - xi) * (1 - et) * (-1 - xi - et - 2 * zt);
		dZt[1] = -0.125 * (1 + xi) * (1 - et) * (-1 + xi - et - 2 * zt);
		dZt[2] = -0.125 * (1 + xi) * (1 + et) * (-1 + xi + et - 2 * zt);
		dZt[3] = -0.125 * (1 - xi) * (1 + et) * (-1 - xi + et - 2 * zt);

		dZt[4] = 0.125 * (1 - xi) * (1 - et) * (-1 - xi - et + 2 * zt);
		dZt[5] = 0.125 * (1 + xi) * (1 - et) * (-1 + xi - et + 2 * zt);
		dZt[6] = 0.125 * (1 + xi) * (1 + et) * (-1 + xi + et + 2 * zt);
		dZt[7] = 0.125 * (1 - xi) * (1 + et) * (-1 - xi + et + 2 * zt);

		dZt[8] = -0.25 * (1 - xi * xi) * (1 - et);
		dZt[10] = -0.25 * (1 - xi * xi) * (1 + et);
		dZt[14] = 0.25 * (1 - xi * xi) * (1 + et);
		dZt[12] = 0.25 * (1 - xi * xi) * (1 - et);

		dZt[9] = -0.25 * (1 - et * et) * (1 + xi);
		dZt[13] = 0.25 * (1 - et * et) * (1 + xi);
		dZt[15] = 0.25 * (1 - et * et) * (1 - xi);
		dZt[11] = -0.25 * (1 - et * et) * (1 - xi);

		dZt[16] = -0.5 * zt * (1 - xi) * (1 - et);
		dZt[17] = -0.5 * zt * (1 + xi) * (1 - et);
		dZt[18] = -0.5 * zt * (1 + xi) * (1 + et);
		dZt[19] = -0.5 * zt * (1 - xi) * (1 + et);
	}

	public static void H20ShapeFunctionCoef(double[][] coord, double[][] coef) {
		double[][] A = new double[20][];
		double[] x = MVO.col(coord, 0);
		double[] y = MVO.col(coord, 1);
		double[] z = MVO.col(coord, 2);
		for (int i = 0; i < 20; i++) {
			A[i] = new double[] {1.0, x[i], y[i], z[i], 
					x[i] * x[i], y[i] * y[i], z[i] * z[i], x[i] * y[i], y[i] * z[i], z[i] * x[i],
					x[i] * x[i] * y[i], x[i] * y[i] * y[i], y[i] * y[i] * z[i], y[i] * z[i] * z[i], 
					x[i] * x[i] * z[i], x[i] * z[i] * z[i], x[i] * y[i] * z[i], 
					x[i] * y[i] * z[i] * x[i], x[i] * y[i] * z[i] * y[i], x[i] * y[i] * z[i] * z[i]};
		}

		IterDSolver solver = new IterDSolver();
		double[] b = new double[20];
		b[0] = 1.0;
		for (int i = 0; i < 19; i++) {
			solver.CGNE(A, b, coef[i], 0);
			b[i] = 0.0;
			b[i + 1] = 1.0;
		}
		solver.CGNE(A, b, coef[19], 0);
	}

	public static void H20ShapeFunction(double[][] xyz, double[] xyz0, double[] N) {
		double[][] coef = new double[20][20];
		H20ShapeFunctionCoef(xyz, coef);
		double x = xyz0[0], y = xyz0[1], z = xyz0[2];
		double[] v = new double[] { 1.0, x, y, z, x * x, y * y, z * z, x * y, y * z, z * x, x * x * y, x * y * y,
				y * y * z, z * z * y, x * x * z, z * z * x, x * y * z, x * x * y * z, x * y * y * z, x * y * z * z };
		MVO.matmul(coef, v, N);
	}

	public static void gradientValue3DH20(double[][] coord, double[] eleU, double[][] gradU) {
		double[][] coef = new double[20][20];
		H20ShapeFunctionCoef(coord, coef);
		double x, y, z;
		double[] dxU = new double[20];
		double[] dyU = new double[20];
		double[] dzU = new double[20];
		for (int i = 0; i < 8; i++) {
			x = coord[i][0];
			y = coord[i][1];
			z = coord[i][2];
			for (int j = 0; j < 20; j++) {
				double[] a = coef[j];
				dxU[j] = a[1] + 2 * a[4] * x + a[7] * y + a[9] * z + 2 * a[10] * x * y + a[11] * y * y
						+ 2 * a[14] * x * z + a[15] * z * z + a[16] * y * z + 2 * a[17] * x * y * z + a[18] * y * y * z
						+ a[19] * y * z * z;
				dyU[j] = a[2] + 2 * a[5] * y + a[7] * x + a[8] * z + a[10] * x * x + 2 * a[11] * x * y
						+ 2 * a[12] * y * z + a[13] * z * z + a[16] * x * z + a[17] * x * x * z + 2 * a[18] * x * y * z
						+ a[19] * x * z * z;
				dzU[j] = a[3] + 2 * a[6] * z + a[8] * y + a[9] * x + a[12] * y * y + 2 * a[13] * y * z + a[14] * x * x
						+ 2 * a[15] * x * z + a[16] * x * y + a[17] * x * x * y + a[18] * x * y * y
						+ 2 * a[19] * x * y * z;
			}
			gradU[i][0] = MVO.dot_product(dxU, eleU);
			gradU[i][1] = MVO.dot_product(dyU, eleU);
			gradU[i][2] = MVO.dot_product(dzU, eleU);
		}
	}

	public static void secondDerivativeValue3DH20(double[][] coord, double[] eleU, double[][] sDU) {
		double[][] coef = new double[20][20];
		H20ShapeFunctionCoef(coord, coef);
		double x, y, z;

		double[] dxxU = new double[20];
		double[] dyyU = new double[20];
		double[] dzzU = new double[20];
		double[] dxyU = new double[20];
		double[] dyzU = new double[20];
		double[] dxzU = new double[20];
		for (int i = 0; i < 20; i++) {
			x = coord[i][0];
			y = coord[i][1];
			z = coord[i][2];
			for (int j = 0; j < 20; j++) {
				double[] a = coef[j];

				dxxU[j] = 2 * a[4] + 2 * a[10] * y + 2 * a[14] * z + 2 * a[17] * y * z;
				dyyU[j] = 2 * a[5] + 2 * a[11] * x + 2 * a[12] * z + 2 * a[18] * x * z;
				dzzU[j] = 2 * a[6] + 2 * a[13] * y + 2 * a[15] * x + 2 * a[19] * x * y;
				dxyU[j] = a[7] + 2 * a[10] * x + 2 * a[11] * y + a[16] * z + 2 * a[17] * x * z + 2 * a[18] * y * z
						+ a[19] * z * z;
				dyzU[j] = a[8] + 2 * a[12] * y + 2 * a[13] * z + a[16] * x + a[17] * x * x + 2 * a[18] * x * y
						+ 2 * a[19] * x * z;
				dxzU[j] = a[9] + 2 * a[14] * x + 2 * a[15] * z + a[16] * y + 2 * a[17] * x * y + a[18] * y * y
						+ 2 * a[19] * y * z;
			}

			sDU[i][0] = MVO.dot_product(dxxU, eleU);
			sDU[i][1] = MVO.dot_product(dyyU, eleU);
			sDU[i][2] = MVO.dot_product(dzzU, eleU);
			sDU[i][3] = MVO.dot_product(dxyU, eleU);
			sDU[i][4] = MVO.dot_product(dyzU, eleU);
			sDU[i][5] = MVO.dot_product(dxzU, eleU);
		}
	}
	
	public static void firstAndSecondDerivativeValue3DH20(double[][] coord, double[] eleU, double[][] fDU, double[][] sDU) {
		double[][] coef = new double[20][20];
		H20ShapeFunctionCoef(coord, coef);
		double x, y, z;
		double[] dxU = new double[20];
		double[] dyU = new double[20];
		double[] dzU = new double[20];
		double[] dxxU = new double[20];
		double[] dyyU = new double[20];
		double[] dzzU = new double[20];
		double[] dxyU = new double[20];
		double[] dyzU = new double[20];
		double[] dxzU = new double[20];
		for (int i = 0; i < 20; i++) {
			x = coord[i][0];
			y = coord[i][1];
			z = coord[i][2];
			for (int j = 0; j < 20; j++) {
				double[] a = coef[j];
				dxU[j] = a[1] + 2 * a[4] * x + a[7] * y + a[9] * z + 2 * a[10] * x * y + a[11] * y * y
						+ 2 * a[14] * x * z + a[15] * z * z + a[16] * y * z + 2 * a[17] * x * y * z + a[18] * y * y * z
						+ a[19] * y * z * z;
				dyU[j] = a[2] + 2 * a[5] * y + a[7] * x + a[8] * z + a[10] * x * x + 2 * a[11] * x * y
						+ 2 * a[12] * y * z + a[13] * z * z + a[16] * x * z + a[17] * x * x * z + 2 * a[18] * x * y * z
						+ a[19] * x * z * z;
				dzU[j] = a[3] + 2 * a[6] * z + a[8] * y + a[9] * x + a[12] * y * y + 2 * a[13] * y * z + a[14] * x * x
						+ 2 * a[15] * x * z + a[16] * x * y + a[17] * x * x * y + a[18] * x * y * y
						+ 2 * a[19] * x * y * z;
				dxxU[j] = 2 * a[4] + 2 * a[10] * y + 2 * a[14] * z + 2 * a[17] * y * z;
				dyyU[j] = 2 * a[5] + 2 * a[11] * x + 2 * a[12] * z + 2 * a[18] * x * z;
				dzzU[j] = 2 * a[6] + 2 * a[13] * y + 2 * a[15] * x + 2 * a[19] * x * y;
				dxyU[j] = a[7] + 2 * a[10] * x + 2 * a[11] * y + a[16] * z + 2 * a[17] * x * z + 2 * a[18] * y * z
						+ a[19] * z * z;
				dyzU[j] = a[8] + 2 * a[12] * y + 2 * a[13] * z + a[16] * x + a[17] * x * x + 2 * a[18] * x * y
						+ 2 * a[19] * x * z;
				dxzU[j] = a[9] + 2 * a[14] * x + 2 * a[15] * z + a[16] * y + 2 * a[17] * x * y + a[18] * y * y
						+ 2 * a[19] * y * z;
			}
			fDU[i][0] = MVO.dot_product(dxU, eleU);
			fDU[i][1] = MVO.dot_product(dyU, eleU);
			fDU[i][2] = MVO.dot_product(dzU, eleU);
			sDU[i][0] = MVO.dot_product(dxxU, eleU);
			sDU[i][1] = MVO.dot_product(dyyU, eleU);
			sDU[i][2] = MVO.dot_product(dzzU, eleU);
			sDU[i][3] = MVO.dot_product(dxyU, eleU);
			sDU[i][4] = MVO.dot_product(dyzU, eleU);
			sDU[i][5] = MVO.dot_product(dxzU, eleU);
		}
	}
	
	public static double[] unitOutNormal3DT3(double[][] xyz) {
		double[] n = new double[3];
		double[] a = {xyz[1][0] - xyz[0][0], xyz[1][1] - xyz[0][1], xyz[1][2] - xyz[0][2]};
		double[] b = {xyz[2][0] - xyz[0][0], xyz[2][1] - xyz[0][1], xyz[2][2] - xyz[0][2]};
		n[0] = a[1] * b[2] - a[2] * b[1];
		n[1] = a[2] * b[0] - a[0] * b[2];
		n[2] = a[0] * b[1] - a[1] * b[0];
		
		double area = Math.sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
		n[0] /= area;
		n[1] /= area;
		n[2] /= area;
		return n;
	}

	/**
	 * 如果点 xyp 在三角形的外接圆内在返回true
	 * 外接圆圆心(xc,yc)，半径r
	 * 注：在外接圆上也算在圆内
	 * @param xyp 点的x坐标
	 * @param x 三角形顶点x坐标
	 * @param y 三角形顶点y坐标
	 * @param circle 三角形外接圆信息, circle[0], circle[1]表示外接圆圆心坐标, circle[2]表示外接圆半径
	 * @return
	 */
	public static boolean CircumCircle(double[] xyp, double[] x, double[] y, double[] circle) {
		double m1,m2,mx1,mx2,my1,my2;
		double dx,dy,rsqr,drsqr;
		double xc, yc, r;
		double x1 = x[0];
		double x2 = x[1];
		double x3 = x[2];
		double y1 = y[0];
		double y2 = y[1];
		double y3 = y[2];
		double EPSILON = 0.000001;
		//检查三点是否重合, 检查x坐标值也可以
		if (Math.abs(y1-y2) < EPSILON && Math.abs(y2-y3) < EPSILON) {
			System.out.println("CircumCircle: Points are coincident.");
			return false;
		}
		//通过计算两个边的中垂线的表达式并求解交点来计算出圆心位置
		//前两种情况为三角形有一条边平行于x轴
		if (Math.abs(y2-y1) < EPSILON) {
			m2 = - (x3 - x2) / (y3 - y2);//中垂线斜率
			mx2 = 0.5 * (x2 + x3);
			my2 = 0.5 * (y2 + y3);
			xc = 0.5 * (x2 + x1);
			yc = m2 * (xc - mx2) + my2;
		}
		else if (Math.abs(y3 - y2) < EPSILON) {
			m1 = - (x2 - x1) / (y2 - y1); //中垂线斜率
			mx1 = 0.5 * (x1 + x2);
			my1 = 0.5 * (y1 + y2);
			xc = 0.5 * (x3 + x2);
			yc = m1 * (xc - mx1) + my1;
		} else {
			m1 = - (x2 - x1) / (y2 - y1);
			m2 = - (x3 - x2) / (y3 - y2);
			mx1 = 0.5 * (x1 + x2);
			mx2 = 0.5 * (x2 + x3);
			my1 = 0.5 * (y1 + y2);
			my2 = 0.5 * (y2 + y3);
			xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
			yc = m1 * (xc - mx1) + my1;
		}
		//计算出半径
		dx = x2 - xc;
		dy = y2 - yc;
		rsqr = dx * dx + dy * dy;
		r = Math.sqrt(rsqr);

		//计算点(xp,yp)到圆心的距离
		dx = xyp[0] - xc;
		dy = xyp[1] - yc;
		drsqr = dx * dx + dy * dy;

		circle[0] = xc;
		circle[1] = yc;
		circle[2] = r;

		return ( drsqr <= rsqr ? true : false );
	}
}
